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We perform some numerical study of the secular triaxial instability of rigidly rotating homogeneous 
fluid bodies in general relativity. In the Newtonian limit, this instability arises at the bifurcation 
point between the Maclaurin and Jacobi sequences. It can be driven in astrophysical systems by 
viscous dissipation. We locate the onset of instability along several constant baryon mass sequences 
of uniformly rotating axisymmetric bodies for compaction parameter M/R — — 0.275. We find 
that general relativity weakens the Jacobi like bar mode instability, but the stabilizing effect is not 
very strong. According to our analysis the critical value of the ratio of the kinetic energy to the 
absolute value of the gravitational potential energy (T/| W|) C rit for compaction parameter as high as 
(N ' 0.275 is only 30% higher than the Newtonian value. The critical value of the eccentricity depends 

very weakly on the degree of relativity and for M/R = 0.275 is only 2% larger than the Newtonian 
value at the onset for the secular bar mode instability. We compare our numerical results with 
recent analytical investigations based on the post-Newtonian expansion. 
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I. INTRODUCTION 
A. The Maclaurin-Jacobi bifurcation point 



CO 
CN 

> 

CN 

O . 

In Newtonian theory a self-gravitating incompressible fluid body rotating at a moderate velocity around a fixed axis 
with respect to some inertial frame takes the shape of a Maclaurin ellipsoid, which is axisymmetric with respect to the 
rotation axis. For a higher rotation rate, namely when the ratio of kinetic to gravitational potential energy T/|W| is 
larger than 0.1375, another figure of equilibrium exists: that of a Jacobi ellipsoid, which is triaxial and rotates around 
its smallest axis JlJ. Actually the Jacobi ellipsoid is a preferred figure of equilibrium, since at fixed mass and angular 
O momentum, it has a lower total energy E = T + W than a Maclaurin ellipsoid, due to its greater moment of inertia / 
with respect to the rotation axis. Indeed, at fixed angular momentum J, the kinetic energy T = J 2 / {21) is a decreasing 
$H . function of I, and for large values of J, this decrease overcomes the effect of the gravitational potential energy W, 
which increases with I. Therefore, provided some mechanism acts for dissipating energy while preserving angular 
momentum (for instance viscosity), a Maclaurin ellipsoid with T/|W| > 0.1375 will break its axial symmetry and 
migrate toward a Jacobi ellipsoid |g, 0, EL This is the secular "bar mode" instability of rigidly rapidly rotating bodies. 
The qualifier secular reflects the necessity of some dissipative mechanism to lower the energy, the instability growth 
rate being controlled by the dissipation time scale 1 . As shown by Christodoulou et al. the Jacobi-like bar mode 
instability appears only if the fluid circulation is not conserved. If on the contrary, the circulation is conserved (as in 
inviscid fluids submitted only to potential forces), but not the angular momentum, it is the Dedekind-like instability 
which develops instead. The famous Chandrasekhar-Friedman-Schutz (CFS) instability (see || for a review) belongs 
to this category. 

The Jacobi-like bar mode instability, applied to neutron stars, is particularly relevant to gravitational wave astro- 
physics. Indeed a Jacobi ellipsoid has a time varying mass quadrupolc moment with respect to any inertial frame, and 
therefore emits gravitational radiation, unlike a Maclaurin spheroid. For a rapidly rotating neutron star, the typical 
frequency of gravitational waves (twice the rotation frequency) falls in the bandwidth of the interferometric detectors 
LIGO and VIRGO currently under construction. 

Neutron stars being highly relativistic objects, the classical critical value T/|W| = 0.1375, established for incom- 
pressible Newtonian bodies, cannot a priori be applied to them. The aim of the present article is thus to investigate 
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For T/|W| > 0.2738, the Maclaurin spheroids are subject to another instability, which is on the contrary dynamical, i.e. it develops 
independently of any dissipative mechanism and on a dynamical time scale (one rotation period). 
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the effect of general relativity on the secular bar mode instability of homogeneous incompressible bodies. We do not 
discuss compressible fluids here. It has been shown that compressibility has little effect on the triaxial instability f| . 

B. Previous studies in the relativistic regime 

Chandrasekhar J?], || has examined the first order post-Newtonian (PN) corrections to the Maclaurin and Jacobi 
ellipsoids, by means of the tensor virial formalism. This work has been revisited recently by Taniguchi |||. However, 
these authors have not computed the location of the Maclaurin-Jacobi bifurcation point at the 1-PN level. This has 
been done only recently by Shapiro & Zane (lO) and Di Girolamo & Vietri jll| . 

On the numerical side, Bonazzola, Frieben & Gourgoulhon [p"2[ |l3| have investigated the secular bar mode instability 
of rigidly rotating compressible stars in general relativity. In the Newtonian limit, they recover the classical result of 
James (see also [Q), namely that, for a polytropic equation of state, the adiabatic index must be larger than 
7 cr it = 2.238 for the bifurcation point to occur before the mass shedding limit (Keplerian frequency). In the relativistic 
regime, they have shown that general relativistic effects stabilize rotating stars against the viscosity driven triaxial 
instability. In particular, they have found that 7 cr it is an increasing function of the stellar compactness, reaching 
7 cr it ~ 2.8 for a typical neutron star compaction parameter. This stabilizing tendency of general relativity has been 
confirmed by the PN study of Shapiro & Zane [jl0| and Di Girolamo & Vietri jll] mentioned above. Note that this 
behavior contrasts with the CFS instability, which is strengthened by general relativity jl(| [l7|] . 

C. The present work 

In this paper, we improve the numerical technique over that used by Bonazzola et al. [|l2| [l3| by introducing 
surface fitted coordinates, which enable us to treat the density discontinuity at the surface of incompressible bodies. 
Indeed the technique used in Refs. JlJ, [l3| did not permit to compute any incompressible model. In particular, it was 
not possible to compare the numerical results in the Newtonian limit with the classical Maclaurin-Jacobi bifurcation 
point. We shall perform such a comparison here. The very good agreement obtained (relative discrepancy ~ 10 -6 ) 
provides very strong support for the method we use for locating the bifurcation point and which is essentially the 
same as that presented in Ref. Jl3| . 

The plan of the paper is as follows. The analytical formulation of the problem, including the approximations we 
introduce, is presented in Sec. |n| Section III then describes the numerical technique we employ, as well as the various 



tests passed by the numerical code. The numerical results are presented in Sec. |Tvj , as well as a detailed comparison 
with the PN studies (lO) and jy]. Finally, Sec. [v| provides some summary of our work. 

II. BASIC ASSUMPTIONS AND EQUATIONS TO BE SOLVED 
A. Helical symmetry 

Let us consider a rotating star that is steadily increasing its rotation rate, e.g. by accretion in a binary system. 
Before the triaxial instability sets in, the spacetime generated by the rotating star can be considered as stationary 
and axisymmetric, which means that there exist two Killing vector fields, k and m, such that k is timelike (at least 
far from the star) and m is spacelike and its orbits are closed curves. 

When the axisymmetry of the star is broken, the stationarity of spacetime is also broken. In Newtonian theory, 
there is no inertial frame in which a rotating triaxial object appears stationary, i.e. does not depend upon the time. 
It can be stationary only in a corotating frame, which is not inertial, so that the stationarity is broken in this sense. 
The stationarity in the corotating frame can be expressed geometrically by stating that the Newtonian spacetime 
possesses a one-parameter symmetry group, whose integral curves are helices. A generator of this symmetry group 
thus has the form 

w 

where t and if are respectively the time and azimuthal coordinates associated with an inertial observer, and f2 is the 
angular velocity with respect to the inertial frame. 

In general relativity, a rotating triaxial system cannot be stationary, even in the corotating frame, as it radiates 
away gravitational waves and therefore loses energy and angular momentum. However, at the very point of the 
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symmetry breaking, no gravitational wave has yet been emitted. For sufficiently small deviations from axisymmetry, 
we may neglect the gravitational radiation. Therefore we shall assume that the spacetime has a helical symmetry, as 
in the Newtonian case. The (suitably normalized) associated symmetry generator £ is then a Killing vector, which 
can be written in the form (^) in weak- field regions (spacelike infinity). Note that spacetimes with helical symmetry 
have been also used for describing binary systems with circular orbits [|l8], |l9[ ^l], |22| . 



B. Rigid rotation 



We model the stellar matter by a perfect fluid, for which the stress-energy tensor takes the form T = (e+p)u<S)u+p g, 
where u is the fluid 4- velocity, e the fluid proper energy density, p the fluid pressure and g the spacetime metric tensor. 

A rigid motion is defined in relativity by the vanishing of the expansion tensor 6 a p := (gj 1 +u a u fJ ")(g j3 l/ -\-upu u )W (vU^ 
of the 4- velocity u. In presence of the Killing vector £, this can be realized by requiring the colinearity of u and £ 
(supposing that the fluid occupies only the region where £ is timelike) : 

u = X£, (2) 

where A is a scalar field related to the norm of £ by the normalization of the 4-velocity A = (— i ■ £)^ 1 ^ 2 . 

For a perfect fluid at zero temperature, the momentum-energy conservation equation V • T = can be recast as 




u • (V A tt) = (3) 



V • (rui) = , (4) 

where n is the proper baryon number density and tt is the momentum I-form tt := hu, h being the fluid specific 
enthalpy: h := (e + p) j \m^n) , where me is some mean baryon mass. In Eq. (|^), VA7T denotes the exterior derivative 
of tt, the so-called vorticity 2-form p3| . For the rigid motion we are considering, the baryon number conservation 
equation (|J) is automatically satisfied, thanks to the colinearity of u with the symmetry generator £. Moreover, the 
equation of motion (|J) can be reduced to a first integral. Indeed Cartan's identity applied to the Lie derivative of the 
I-form tt along the vector field £ leads to 

£ t ix = £ ■ (V A tt) + V(£ ■ tt) = , (5) 

where the second equality holds thanks to the helical symmetry: £ grr = 0. Inserting relation (^) into the equation of 
fluid motion ||) shows that the first term in Eq. (Q) vanishes identically, so that one gets the first integral of motion 
11 

£ ■ tv = —X^ 1 h — const. (6) 



In the axisymmetric and stationary case, where £ is a linear combination of the two Killing vectors [Eq. (14) below], 



one recovers the classical expression 25 



The first integral (|6|) can be re-expressed in terms of the 3+1 formalism of general relativity. Let us introduce the 
spacetime foliation by the t = const hypersurfaces St and the associated future-directed unit vector field n everywhere 
normal to Ej. n is the 4-velocity of the so-called Eulerian observer (also called ZAMO or locally non-rotating observer). 
We have the following orthogonal split of the fluid 4-velocity: 

u = r(n + U) with n U = , (7) 

with the Lorentz factor 

r= n u = (l U U) 1/2 . (8) 

The spacelike vector U is the fluid 3-velocity as measured by the Eulerian observer and the second equality in the 
above equation results from the normalization of the 4-velocity u. Similarly, we have an orthogonal split of the helical 
Killing vector: 



£ = Nn + B with n B = , 



(9) 
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where TV is the lapse function, governing the proper time evolution between two neighboring hypersurfaces £4. From 
relation (2) and the two orthogonal decompositions (0) and (||), we get A = T/N, so that the (logarithm of the) first 
integral (6) can be written 

H + u-lnT = const. (10) 

with 

H:=\nh and ^:=lniV. (11) 

In the non-relativistic limit, H tends toward the classical specific enthalpy (excluding the rest-mass energy) of the fluid, 
whereas v tends toward the Newtonian gravitational potential. Therefore, in the Newtonian limit, where lnT — > U 2 /2 
and U — > SI x r, Eq. (|l0|) reduces to the classical first integral of motion 

H + v - i(n x r) 2 = const. (12) 
C. Einstein equations 

When the rigidly rotating star is still stationary and axisymmetric, it is well known that a coordinate system 
x^ = (t, r, 9, ip) can be chosen so that the metric takes the Papapetrou form (see e.g. (^(| and references therein) 

dx»dx v = -N 2 dt 2 + B 2 r 2 sin 2 6(d<p - N v dt) 2 + A 2 (dr 2 + r 2 d(9 2 ) , (13) 

where N, N v , A and B are four fu nctio ns of (r, 9). The coordinate vectors d/dt and 8/ dip are then the two Killing 



vectors k and m mentioned in Sec. II A . The helical Killing vector t is expressible as 

^=k + f2m. (14) 



With the form (13), the Einstein equations reduce to a set of four coupled elliptic equations (see e.g. [|27| for the 
precise form). 

When the triaxial instability sets in, we shall consider that the metric is a perturbation of (IKS)), which we write as 

g^y dx^dx v = -N 2 dt 2 + B 2 r 2 sin 2 6(d<p - N v dt) 2 + A 2 (dr - N r dt) 2 + r 2 A 2 (d9 - N e dt) 2 , (15) 

where N, N r , N e , A and B are six functions of (r, 9, ip'), with 

ip' := <p - Sit . (16) 

Note that 8/dt and 8 /dip are no longer Killing vectors, only I remains Killing. The coordinate system x M = (t, r, 9, ip') 
is adapted to the Killing vector £ and t is an ignorable coordinate in this system. We call x^ the corotating coordinate 
system, and x^ the nonrotating one. N is the lapse function already introduced in Eq. (^|). N l = (N r , N e , N v ) is 
(minus 2 ) the shift vector of the nonrotating coordinate system. 

In the triaxial case, there is no equivalent of the Papapetrou theorem |2^, which, in absence of convective 

motions (in the meridional planes), allows one to set to zero all the off-diagonal components of the metric tensor 
of axisymmetric and stationary spacetimes, except for g tlp . So, in principle, all the metric components should be 
non- vanishing in Eq. (^|). However, we retained only g tr and gte as the extra non zero components with respect 
to the axisymmetric case. We did so as an approximation in order to simplify the writing of Einstein equations. 
This approximation can be justified by the following remarks: (i) the metric element (|l^) is exact in the stationary 
axisymmetric case, (ii) it encompasses the first order PN metric. This means that the neglected terms are of second 
order PN and moreover vanish for stationary axisymmetric rotating stars. We consider these terms to be negligible 
in our study of the verge of the non-axisymmetric instability. 



2 In numerical relativity, the current convention is to call shift vector the 3- vector /3 l = 



-N\ 
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The Einstein equation leads to the following set of partial differential equations: 



Au = 4TrA 2 (E + 3p+{E+p)U t lf l )+A 2 K lJ K lj -V l vV\v + l3) (17) 
A AT 4 + -V*V,-iV j = -16ttNA 2 {E +p)U l + NB- 2 K lj W 3 (6(3 - v) (18) 
A 2 [(NB - l)rsin0] = 167riV A 2 Bp r sin 9 (19) 



A 2 C = SttA 2 [p+{E + p)U t U % ] + -A 2 K l jK i ° -VivV v , (20) 
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where the following notations have been introduced [see also Eq. (|TT|)1 

C := ln(AZV) , /3:=lnB , (21) 

and V, denotes the covariant derivative with respect to the flat 3-metric = diag(l, r 2 , r 2 sin 2 6), A = V^V* the 
corresponding Laplacian and A 2 is the Laplacian in the 2-dimensional space spanned by (r, 6) : 

A -— + -— + —— (22) 
dr 2 r dr r 2 d6 2 

In Eqs. (p7|)-(|20|), U l denotes the spatial components of the fluid 3- velocity U [cf. Eq. ([?])]. They can be expressed as 

N r N e 1 

U r = , U B = , U v = — (Q-N*), (23) 

N N N y ' 

and r is deduced from U by means of Eq. (|J) : 

r = (1 - A 2 [{U r f + r 2 {U 6 ) 2 ] - B 2 r 2 sin 2 6(U^) 2 y 1/2 . (24) 
We have also introduced the energy density E as measured by the Eulerian observer: 

E = n • T- n = T 2 (e+p) -p . (25) 



Another quantity not yet defined and which appears in Eqs. (|17|)-(|20|) is the extrinsic curvature tensor K of the 
hypersurface Sj. 

D. Equation of state 

The above equations must be supplemented by an equation of state (EOS), i.e. a relation between H , which appears 
in the first integral ([l(]), and e and p which appear in the source terms of the Einstein equations jl7|)-(|20|), via Eq. (pH). 
For the incompressible matter we are considering, we choose this EOS to be 

e(H) = po (26) 
p(H) = p (exp(ff) - 1) (27) 
n(H) - Po /m B , (28) 

where po is a constant, representing the constant proper energy density of the fluid. 

III. NUMERICAL CODE AND ITS TESTS 
A. Numerical technique 

The basic idea is to solve the equations presented in Sec. |l| by means of an iterative scheme to get an axisymmetric 
stationary solution, and then to introduce some triaxial perturbation and resume the iterative scheme. If the pertur- 
bation is damped (resp. grows) as long as the iteration proceeds, the equilibrium configuration will be declared stable 
(resp. unstable). This procedure has been used already in the works |l^, |^|. We shall prove here, by comparison 
with the analytical result, that it correctly locates the secular instability point along the Maclaurin sequence. 
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FIG. 1: Evolution during the iterative procedure of the convergence indicator SH/H (long-dashed line), the triaxial perturba- 
tion q (short-dashed line) and its growth rate q (solid line). The left figure corresponds to a stable configuration with respect 
to nonaxisymmetric perturbation, the right one an unstable one. SH/H and q are depicted in logarithmic units, while q is 
multiplied by 10 for better readability. The discontinuity at step 10 is due to the switch of the rotation and that at step 25 to 
the switch of the triaxial perturbation. 



We solve the system of non-linear elliptic equations (J17|)-(j20|) by means of the multi-domain spectral method 
presented in Ref. |3l). The nice feature of this method is that it introduces surface-fitted coordinates (£,0',<p'), so 
that the density discontinuity at the stellar surface is exactly located at the boundary between two domains. In this 
way, all the fields are C°° functions in each domain. This avoids any spurious oscillations (Gibbs phenomenon) and 
results in a very high precision. As discussed in Sec. |[ Cj , this technique constitutes the major improvement with 
respect to the numerical method used in Ref. Jll|. Another difference with Ref. [jl3| is that we employ the technique 
presented in Ref. to solve the vector elliptic equation ( |l8| ) for the shift vector. In particular, we solve for the 
Cartesian components of the shift vector, whereas |l3| solved for the spherical components. 

A weak point of the surface-fitted coordinate technique of Ref. |3l[ is that the transformation from the computational 
coordinates (£, 6' , ip') to the physical ones (r, 8, p) becomes singular when the ratio of polar to equatorial radius r p /r eq 
is lower than ~ 1/3. This prevents us from computing very flattened configurations, but fortunately this causes no 
trouble in reaching the Maclaurin-Jacobi bifurcation point, which is at r p /r cq ~ 0.58. 

The iterative procedure is as follows. First one must set a value of the central enthalpy H c and the rotation 
frequency f2, in order to pick out a unique rotating axisymmetric configuration. The iteration is then started from 
very crude values: flat spacetime and spherical stellar shape. We solve the Einstein equations (|17|)-(po|) with the 
spherical energy density and pressure distribution in their right-hand side. We then plug the obtained value of v in 
the first integral of motion |h]), along with H c , to get the enthalpy field H. The zero of this field defines the new 
stellar surface to which we adapt the computational coordinates (£, 6', ip') (s ee pl| for details). From H we compute 
new values of the energy density and pressure according to the EOS (pq)- (p7|)7These values are then put on the 
right-hand side of the Einstein equations ([l7])-(^0]) and a new iteration begins. Usually the rotation velocity is set 
to zero for the ten first steps. It is then switched on, either integrally or (for very rapidly rotating configurations) 
gradually. The convergence of the procedure is monitored by computing the relative difference SH/H between the 
enthalpy fields at two successive steps (long-dashed line in Fig. Q). The iterative procedure is stopped when SH/H 
goes below a certain threshold, typically 10~ 7 , or for high precision computations 10~ 12 . 

After this convergence has been achieved, we switch on a triaxial m = 2 perturbation by modifying the metric 
potential v according to 

^^j/(l + esin 2 6»cos2( / 3) , (29) 

where e is a small constant, typically e ~ 10~ 5 to 10~ 3 . Then we continue the iteration procedure as described above, 
without any further modification of the equations. At each step, we evaluate the quantity 



q := max p% + pf , (30) 
where v% and vj, are the m — 2 coefficients of the Fourier expansion of the tp part of v. 



v(tp) = vq + z>2 cos(2<p) + i>3 sin(2iy9) + £4 cos(4cp) + £5 sin(4y>) + 



(31) 
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FIG. 2: Isocontours of the enthalpy H in the equatorial plane at steps 20, 150, 180 and 200 for the unstable model presented 
in Fig. [j]. Dashed lines denote negative values of H, which correspond to the exterior of the star. The thick solid line denotes 
the stellar surface. 



and the max in Eq. ([50|) is taken over all the r and 9 coefficients, q is used to monitor the evolution of the triaxial 
perturbation: if q — > as the iteration proceeds, we conclude that the perturbation decays and the axisymmetric 
configuration is a stable one. In the vicinity of the marginally stable configuration, the decay or growth of the 
perturbation turns out to be pretty small. In order to facilitate the diagnostic, we monitor instead the relative growth 
rate of the m = 2 part, defined as 

Q ^previous step /or*\ 

q := ^ z . (32) 

^previous step 

After some transitory regime, q turned out to be constant. If q > (resp. q < 0) we conclude that the configuration 
is unstable (resp. stable). The behaviors of SH/H, q and q in two typical computations are shown in Fig. [j], whereas 
Fig. H depicts the development of the triaxial instability. 

Note that the above method does not require to specify some value of viscosity. Whatever this value, the effect of 
viscosity is simulated by the rigid rotation profile that we impose at each step in the iterative procedure via Eq. (|2^). 
If we consider the iteration as mimicking some time evolution, this means that the time elapsed between two successive 
steps has been long enough for the actual viscosity to have rigidified the fluid flow. Consequently, a quantity that we 
cannot get by our method is the instability time scale. As recalled in the Introduction, this time scale depends upon 
the actual value of viscosity, but not the instability itself. 

The numerical code implementing the above procedure has been constructed upon the C++ library Lorene [j33| . 
Numerical computations have been performed on SGI Origin200 as well as Linux PC workstations. Three domains 
have been used, the innermost one corresponding to the interior of the star, and the outermost one extending to 
spacelike infinity by means of a suitable compactification (see for details). The number of spectral coefficients 
used in each domain is N r x Ng x N v — 33 x 17 x 4. The corresponding memory requirement is 40 MB and a typical 
CPU time (e.g. corresponding to the first 60 steps of Fig. [I], which are sufficient to conclude about the stability of 
the configuration) is 3 min on an Intel Pentium IV 1.5 GHz processor. 



B. Test of the numerical code 



Numerous tests have been performed to assess the validity of the method and the accuracy of the numerical code. We 
present here successively tests for axisymmetric configurations in general relativity and tests about the determination 
of the triaxial insta bility point in the Newtonian limit. Tests regarding the triaxial instability in the rclativistic case 
are deferred to Sec. IV B, where we present a detailed comparison with analytical (post-Newtonian) results. 



1. Tests in the axisymmetric regime 



The multi-domain spectral technique with surface-fitted coordinates has already been tested in the Newtonian 
regime, giving a rapidly rotating Maclaurin ellipsoid with a relative error of the order of 10 -12 It has been 

also shown in Ref. |5lJ that the error is evanescent, i.e. that it decays exponentially with the number of spectral 
coefficients, which is typical of spectral methods. 
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TABLE I: Comparison with the numerical results of Ansorg et al. [|38|] for a highly relativistic stationary axisymmetric model 
Pc/po = 1 (H c = hi2) and r p /r cq = 0.7, with compaction parameter M/R c i IC ~ 0.39. Meaning of the symbols [in geometrized 
units (G — c — 1)] are as follows: dimensionless angular velocity Q := Q,/ p]J 2 , gravitational mass M := Mp]/ 2 , baryon mass 
Mb := MbP , circumferential equatorial radius ii c i rc := RdicPo , angular momentum J := Jpo- Z p , Z BCl , and Z° q are the 
redshifts in the polar direction, in the equatorial plane along the direction of motion and opposite to that direction respectively. 





our value 


relative diff. 


n 

~M 

M 

-^circ 

J 

z P 

7 { 

^cq 
yb 
^cq 


1.4116964 

0.1358182 

0.1863648 

0.34548954 

0.01405836 

I. 7073708 
-0.1625613 

II. 3539974 


0.0009 % 
0.015 % 
0.014 % 
0.0039 % 
0.0017 % 
0.001 % 
0.017 % 
0.00073 % 



The accuracy of the computed relativistic axisymmetric models is estimated using two general relativistic virial 
identities GRV2 |34|, [35] and GRV3 |3(| . These two virial error indicators are integral identities which must be satisfied 
by any solution of the Einstein equations and which are not imposed during the numerical procedure (see Ref. 37 for 
the computation of GRV2 and GRV3). GRV3 is a generalization to general relativity of the classical virial theorem. 
We have obtained values of GRV2 and GRV3 of the order of 10~ 12 in Newtonian regime and 10 -8 — 10~ 4 for high 
rotation and high compaction parameter. The virial errors are at least one order of magnitude better than obtained 
in the previous code used for calculating rapidly rotating strange stars |}9[ [l0[ [ll], |^|, ^3], [44|, |^|, |iq| . 

Another test in the strongly relativistic regime is the comparison with the highly accurate recent code of Ansorg 



et al. 1 38 , also based on spectral methods but using very different coordinates and resolution scheme. The relative 
differences defined as diff = |(AKM - GG)/AKM| (where AKM and GG are the values obtained by Ansorg et al. ||| 
- their Table 1 - and us respectively) , are presented in Table | for an incompressible rapidly rotating body with a 
compaction parameter M/R C - 1IC — 0.393, where i? c ; rc is the circumferential radius of the star. Note that this model is 
much more relativistic than any realistic neutron star, so that we are really testing the strong field behavior of our 
code. Table | shows that the agreement with Ansorg et al. is very good, of the order of 10~ 4 or better. For this 
model, the GRV2 and GRV3 errors from our code are 2 10 -5 and 3 10 -6 respectively. 



2. Tests about the code capability to find the triaxial instability 

Widely used indicators for the onset of instability are the eccentricity e and the ratio of the kinetic energy to the 
absolute value of the gravitational potential energy T/lW^I, defined respectively as 

e 2 = 1 - (r p /r oq ) 2 , (33) 

T/\W\ = nJ ^ 2 (34) 

71 1 fiJ/2 + M B -M' 1 ' 

where r p and r eq are the polar and equatorial coordinate radius, J the total angular momentum and Mb the baryon 
mass. Note that is gauge invariant for uniformly rotating incompressible objects 3 , while eccentricity is coor- 

dinate dependent. In the Newtonian regime, we have employ the code to locate the bifurcation point between the 
Maclaurin sequence and the Jacobi one, and found it to be at 

(T/\W\) crit = 0.137526 and e crit = 0.812667 . (35) 

Both quantities differ from the exact value (T , /|W / |) cr j ti N ewt and e cr i tj Ncwt by only 10~ 6 . 



3 For compressible bodies, T/|W| can be defined according to Eqs. (20)- (22) of Ref. |47| , for which M-g at the denominator of Eq. jj^ ) 
must be replaced by some "proper" internal energy M p ; in general relativity this last quantity has less physical meaning than the total 
baryon mass Mj. 
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FIG. 3: Square of the angular velocity (in units of nGpa) as a function of the eccentricity [defined by Eq. (p3|)] for axisymmetric 
equilibrium sequences of constant rest mass. Each sequence is labeled by the compaction parameter M s /R s of its nonrotating 
member. The dotted line denotes the Maclaurin sequence {M B /R S = 0). The thick line connects the secular instability points. 
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FIG. 4: Square of the angular velocity versus the ratio of the kinetic energy to the absolute value of the gravitational potential 
energy [defined by Eq. (|34j)] for the same sequences as in Fig. |§| The thick line connects the secular instability points. 
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FIG. 5: Comparison between the present relativistic computations of the bar mode instability points (thick solid line with 
circles) and that of two different PN calculations performed by 0] (thick solid line with triangles) and by fll| (thick solid line). 
The Newtonian Maclaurin-Jacobi bifurcation point is shown as a filled circle. Thin solid lines and long-dashed lines corresponds 
respectively to fully relativistic calculations and to the ellipsoidal PN ones of (l^] for axisymmetric equilibrium sequences with 
compaction parameter M s /R s = 0.075; 0.05; 0.025 from top to bottom. The dotted line represents the Newtonian Maclaurin 
sequence (M B /R B = 0). 



IV. CALCULATIONS AND RESULTS 
A. Rotating axisymmetric equilibrium configurations 

We present in this section our resu lts about the relativistic analogs of the axisymmetric Newtonian Maclaurin 



ellipsoids. As mentioned in Sec. Ill A, the multi-domain spectral technique with surface-fitted coordinates that we 
employ has some trouble when r p /r oq < 1/3. To compute very rapidly rotating configurations, well beyond the 
Jacobi-like bar mode instability (which is located at r p /r cq ~ 0.58 in the Newtonian reg ime), we employ a second 
numerical code for axisymmetric stationary relativistic stars, namely that of Stergioulas ||48f (see |fi9f for a description). 
More precisely, the multi-domain spectral code has been used for calculating rotating stars with the eccentricity e (or 
T/IM^I) lower than 0.92 — 0.87 (0.23-0.2) for compaction parameters ranging from to 0.25, while Stergioulas' code 
has been used for higher values. In the overlapping region around e ~ 0.9, models computed by both codes agree very 
well in all computed properties. 

We have constructed constant baryon mass sequences of uniformly rotating axisymmetric homogeneous fluid bodies. 
Each of the sequences is parametrized by the compaction parameter M s /R s , where M s and R s are the gravitational 
mass and circumferential radius of the spherical (non-rotating) member of the sequence. To show the role of relativistic 
effects we plot in Fig. || and Fig. || the square of the angular velocity as a function of the eccentricity e and T/|W| 
respectively for several axisymmetric equilibrium sequences (thin lines) with compaction parameters M s /R s ranging 
from (Newtonian Maclaurin sequence) to 0.25. We find that in general relativity the rotational velocity along 
the equilibrium sequence is not much different from that one derived in the Newtonian limit, although a star of a 
given eccentricity has a little bit higher value of the square of the angular velocity (the biggest difference is near the 
maximum, and for M s /R s = 0.25, it is by 35% higher than the Newtonian one). The location of the maximum of 
n 2 /(7rp ) weakly depends on the compaction parameter and is only a slightly shifted to lower values of eccentricity 
and higher values of T/|W| for highly relativistic stars. 
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FIG. 6: Square of the angular velocity as a function of the ratio of the kinetic energy to the absolute value of the gravitational 
potential energy for the same equilibrium sequences as in Fig. [| (all symbols and lines have the same meaning as in Fig. ^). The 
thick lines with circles and filled triangles correspond to the secular instability points obtained by us and by [[To| respectively. 

In Figs. |5| and ^ we compare our fully relativistic sequences with PN equilibrium sequences in the ellipsoidal 
approximation |Q, for compaction parameters M s /R s — 0.025, 0.05 and 0.075. Let us recall that, contrary to the 
Newtonian case, in general relativity the figure of equilibrium of rotating incompressible bodies is not an ellipsoid 
(with respect to the employed coordinates). However as it can be seen for low compactness the ellipsoidal PN models 
and fully relativistic models are in close agreement up to the value of eccentricity ~ 0.9 and r/|W| ~ 0.2. The 
discrepancy increases as eccentricity increases. For M s /R s = 0.075 and e = 0.9 the value of f2 2 /(ftpo) is overestimated 
by 3% using the ellipsoidal approximation. 

For high compaction parameter (M s /R s > 0.1) comparison between our fully relativistic calculations (Figs. ||and^) 
and PN ones (Fig. 3 in [ fl0| ) shows that the ellipsoidal approximation does not provide a good description of equilibrium 
axisymmetric configurations. We find a more pronounced increase of the angular velocity with compactness at any 
given value of eccentricity. The relative differences are 10% and 20% for Af s /i? s = 0.15 and 0.25 respectively and 
e < 0.8 (there are much larger for higher value of e). Moreover according to |l0| the location of the maximum of 
f2 2 /(7rpo) is shifted to higher values of e with respect to Newtonian one, while we find the opposite tendency, in 
agreement with JTl[ ] . 

B. Secular bar mode instability 

1. Our results 

Some important quantities at the viscosity driven instability points for different values of compaction parameter 
are reported in Table [fl| We are using two kinds of compaction parameters: M s / R s already defined and the proper 
compaction parameter M/R c i rc , where i? c i rc is the circumferential radius, i.e. the length of the equator (as given by 
the metric) divided by 2tt of the actual configuration (unlike R s which is the circumferential radius of the non rotating 
star having the same baryon mass). 

Our results are shown as thick lines in Figs. [| || and|7|. According to our analysis the critical value of the eccentricity 
very weakly depends on the compaction parameter and for M s /i? s = 0.275 is by only 2% larger than Newtonian value 
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TABLE II: Jacobi-like instability points along relativistic sequences. The symbols are as follows (all values are for 
axisymmetric configurations at the instability point): H c is the central enthalpy; fi M is the proper compaction parameter, 

where -Rcirc is circumferential radius; is the compaction parameter of nonrotating spherical configurations with the same rest 
mass as the marginally stable configuration; N c is central lapse; f2 2 /(7rpo) is a square of the angular velocity in the unit ivGpo; 
e cr it and (T/|W|) C rit are eccentricity and the ratio of the kinetic energy to the absolute value of the gravitational potential 
energy at the instability points; GRV2 and GRV3 are virial errors. 
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0.01 
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-3e-07 
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5e-05 


0.04 
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0.15097 
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0.65885 


0.49528 
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0.2000 
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0.15 
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0.53192 


0.82795 
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2e-05 


3e-05 


0.16 


0.258 


0.2272 


0.55757 


0.54057 
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0.17859 


le-05 


-6e-05 


0.18 


0.275 


0.2430 


0.52322 


0.55748 


0.82838 


0.18220 


-3e-05 


-2e-4 



of the onset of the secular bar mode instability. The critical value of the ratio of the kinetic energy to the absolute 
value of the gravitational potential energy (T/\ W\) cr i t for compaction parameter as high as 0.275 is by 30% (Table ||) 
higher than the Newtonian value. The dependence of (T/|W|) cr i t on the compactness can be very well approximated 
by the function 

frri\tim w;ni7ii i J 0.148x(x + 1) for x:=M/R circ , 

(T/\W\U t ^(T/\W\) crm+ 0J2Mx+1) for x:=Ms/Rs (36) 



2. Comparison with PN calculations 



The comparison between our relativistic calculations of the viscosity driven instability and the corresponding two 
different PN calculations derived by 10 and by Jll| is shown in Fig 0, || and |[ Shapiro and Zane jl(| employ 
an energy variational principle to determine equilibrium shape and stability of homogeneous triaxial ellipsoids. The 
method used by them is valid for arbitrary rotation rate, but only for constant density bodies. Di Girolamo and 
Vietri |ll[ determined the value of the eccentricity and fl 2 /(irpo) at the instability onset point using Landau's theor 
of second-order phase transitions. This method is the extension to PN regime of that used by Bertin and Radicati 
for the Newtonian treatment of bar mode instability and valid for any equation of state. Considering the dependence 
of fl 2 1 (irpo) with respect to the eccentricity (Fig. ||) we found good agreement between our results and the PN ones of 
jllj . As can be seen, both calculations show much weaker influence of relativistic effects on location of the instability 
onset point than Shapiro & Zane PN calculations. 

Figure || shows the dependence of the critical eccentricity and the critical value of T/|VF| on the compaction 
parameter M s /R s . The Authors of [jllj use a proper "conformal compaction parameter" M c /R c := Anpo/SR 2 , 
where R c = (aiazCLs) 1 ^ , Qi, 02, 13 being the ellipsoid semiaxes. In order to make a comparison between our 
relativistic calculations and those of |ll| we find the M a /R s corresponding to their M c /R c using the formula M c /R c = 
0.25M s /i? s (l - M s /R s + (1 - 2M S /R S )) 2 (see f(j, expression [4], p. 422). This formula is valid in the spherical limit, 
but the difference between the proper "conformal compaction parameter" and the "spherical conformal compaction 
parameter" is at most 3 %. 

According to our analysis the critical value of the eccentricity depends very weakly on the compaction parameter 
(solid line); this is in agreement with the PN calculations by Di Girolamo & Vietri The relative differences 

between our and their PN calculations are less than 1%. A much stronger weakening of the bar mode instability by 
general relativity is suggested by the PN study of Shapiro & Zane jn^, who find that e cr it (dashed line in the left 
panel of Fig. ||) could be as large as 0.94 at M s /R s = 0.25. This discrepancy may be ascribed mainly to the ellipsoidal 
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FIG. 8: Eccentricity (left) and ratio of the kinetic energy to the absolute value of the gravitational potential energy (right) at 
the onset of the secular bar mode instability, versus the compaction parameter M a /Rs. Our results are denoted by the solid 
line. The dashed line and the dash-dotted line correspond to the PN calculations by mXX and fill] respectively. 



approximation for the deformation and equilibrium shape. We refer to the article |ll| for a detailed explanation of 
discrepancies between the two PN calculations of instability points. 

The right panel of Fig. || presents the comparison between |l(| and our calculation of the critical T/|W| as a 
function of the compaction parameter M s /R s . We found that for the compaction parameter as high as 0.05; 0.15; 
0.25 (T/|W|) cr it is only by 5%; 16%; 28% higher than Newtonian value, while according to Ref. [jl0|, the increase is 
40% ; 79% ; 88% respectively. 

According to our study relativistic effects weaken the Jacobi-like bar mode instability (solid line in Figs. |[ and 
|8|), but the stabilizing effect is not very strong. This is in agreement with the PN calculations pT| . 



V. SUMMARY 



Triaxial instabilities of rotating compact stars can play an important role as emission mechanisms of gravitational 
waves in the frequency range of the forthcoming interferometric detectors. A rapidly rotating neutron star can 
spontaneously break its axial symmetry if the ratio of the rotational kinetic energy to the absolute value of the 
gravitational potential energy T/| W| exceeds some critical value. We have investigated the effects of general relativity 
upon the nonaxisymmetric viscosity-driven bar mode instability of incompressible, uniformly rotating stars. This is 
the relativistic analog of the Newtonian Maclaurin-Jacobi bifurcation point. 
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Our method of finding the instability point is similar to that used in Ref. []13| for compressible fluid stars. The main 
improvement with respect to this work regards the numerical technique. We have indeed introduced surface-fitted 
coordinates, which enable us to treat the strongly discontinuous density profile at the surface of incompressible bodies. 
This avoids any Gibbs-like phenomenon and results in a very high precision, as demonstrated by comparison with the 
analytical result for the Newtonian Maclaurin-Jacobi bifurcation point, that our code has retrieved with a relative 
error of 10~ 6 . 

According to our results, general relativity weaken the Jacobi-like bar mode instability: the values of T/W, ec- 
centricity and fl 2 /(irpo) increase at the onset of instability above the Newtonian values. This general tendency is 
in agreement with PN analytical results |l^, [ll[] for rigidly rotating incompressible bodies and with the numerical 
calculations of |u| for relativistic polytropes. However we found that the stabilizing effect of general relativity is 
much weaker than that obtained in the PN treatment of fTH] . The critical value of the ratio of the kinetic energy to 
the absolute value of the gravitational potential energy (T/IPTDcrit f° r a compaction parameter as high as 0.275 is 
only 30% larger than the Newtonian value, whereas it has been found 90% larger by (l(J. 

According to our analysis the critical value of the eccentricity very weakly depends on the compaction parameter 
and for a compaction parameter as high as 0.275 is only 2% (16% according to Jl_0| ) larger than the Newtonian value 
of the onset of the secular bar mode instability. Regarding the dependence of 2 /(7rpo)(e cr it) and e cr it with respect to 
compactness, we found very good agreement between our result and the recent PN ones of ||ll| , the relative differences 
being lower than 1%. 
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